stanEx4
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # not see parameters str1..., str2,... using regex as exc='[str1,str2,...]'
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(4,1),mar=c(1,5,1,1))
drw[,1,] |> plot(type='l',xlab='',ylab=pr)
drw[,2,] |> plot(type='l',xlab='',ylab=pr)
drw[,3,] |> plot(type='l',xlab='',ylab=pr)
drw[,4,] |> plot(type='l',xlab='',ylab=pr)
par(mar=c(3,5,3,3))
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)
mcmc=goMCMC(mdl,data,smp,wrm)
mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
}explanatory variable obs.x have noise
xN(x0.sx),yN(b0+b1*x0,s)
ex8-0.stan
normal regression without explanatory variable’s noise
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}ex9.stan
normal regression with explanatory variable’s noise
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x10;
}
parameters{
real b0;
real b1;
real<lower=0> s;
real<lower=0> sx;
vector[N] x0;
}
model{
for(i in 1:N){
x[i]~normal(x0[i],sx);
y[i]~normal(b0+b1*x0[i],s);
}
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x0[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] x1;
vector[N1] y1;
for(i in 1:N1){
x1[i]=normal_rng(x10[i],sx);
m1[i]=b0+b1*x10[i];
y1[i]=normal_rng(m1[i],s);
}
}n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)
n1=10
#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -9805.7
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 26 -38.3056 0.000609277 0.000161797 1 1 59
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -38.31
## b0 8.48
## b1 1.71
## s 4.12
## m0[1] 37.20
## m0[2] 19.52
## m0[3] 27.93
## m0[4] 7.43
## m0[5] 25.25
## m0[6] 42.66
## m0[7] 38.88
## m0[8] 22.16
## m0[9] 43.71
## m0[10] 23.40
## m0[11] 4.86
## m0[12] 38.17
## m0[13] 46.06
## m0[14] 32.24
## m0[15] 36.19
## m0[16] 20.36
## m0[17] 25.05
## m0[18] 41.02
## m0[19] 14.27
## m0[20] 31.03
## y0[1] 35.78
## y0[2] 27.80
## y0[3] 26.60
## y0[4] 8.80
## y0[5] 23.02
## y0[6] 39.06
## y0[7] 40.91
## y0[8] 20.69
## y0[9] 42.88
## y0[10] 23.26
## y0[11] 4.09
## y0[12] 44.86
## y0[13] 56.49
## y0[14] 31.82
## y0[15] 35.52
## y0[16] 23.46
## y0[17] 25.41
## y0[18] 41.26
## y0[19] 16.71
## y0[20] 34.24
## m1[1] 4.86
## m1[2] 9.44
## m1[3] 14.02
## m1[4] 18.59
## m1[5] 23.17
## m1[6] 27.75
## m1[7] 32.33
## m1[8] 36.91
## m1[9] 41.48
## m1[10] 46.06
## y1[1] 1.90
## y1[2] 9.00
## y1[3] 18.24
## y1[4] 19.46
## y1[5] 20.95
## y1[6] 26.51
## y1[7] 34.63
## y1[8] 40.97
## y1[9] 42.02
## y1[10] 43.58
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -38.58 -38.23 1.38 1.16 -41.38 -37.10 1.01 421 392
## b0 8.37 8.36 2.29 2.04 4.47 12.03 1.02 272 264
## b1 1.71 1.72 0.17 0.16 1.45 1.98 1.01 346 426
## s 4.66 4.56 0.85 0.79 3.54 6.21 1.00 1167 1252
## m0[1] 37.23 37.24 1.36 1.37 34.97 39.37 1.00 2056 1473
## m0[2] 19.46 19.48 1.46 1.36 17.07 21.76 1.01 301 283
## m0[3] 27.91 27.88 1.12 1.09 26.09 29.76 1.01 610 889
## m0[4] 7.31 7.32 2.38 2.10 3.26 11.15 1.02 272 265
## m0[5] 25.21 25.19 1.18 1.11 23.26 27.16 1.01 427 528
## m0[6] 42.70 42.70 1.71 1.63 39.81 45.43 1.00 1270 1119
## m0[7] 38.91 38.92 1.45 1.47 36.46 41.19 1.00 1811 1410
## m0[8] 22.11 22.11 1.31 1.22 19.93 24.21 1.01 337 373
## m0[9] 43.77 43.76 1.79 1.70 40.74 46.66 1.00 1157 965
## m0[10] 23.35 23.34 1.25 1.18 21.27 25.36 1.01 364 417
## m0[11] 4.73 4.77 2.61 2.33 0.40 8.92 1.02 273 268
## m0[12] 38.20 38.21 1.41 1.43 35.82 40.40 1.00 1929 1377
## m0[13] 46.13 46.11 1.98 1.86 42.78 49.26 1.00 967 923
## m0[14] 32.24 32.24 1.15 1.14 30.35 34.10 1.00 1420 1308
## m0[15] 36.21 36.23 1.30 1.32 34.04 38.28 1.00 2123 1361
## m0[16] 20.30 20.32 1.41 1.32 17.98 22.56 1.01 310 322
## m0[17] 25.02 25.00 1.18 1.12 23.05 26.97 1.01 418 520
## m0[18] 41.06 41.06 1.60 1.55 38.38 43.55 1.00 1477 1331
## m0[19] 14.18 14.19 1.83 1.63 11.08 17.03 1.02 277 221
## m0[20] 31.02 31.01 1.13 1.12 29.20 32.86 1.00 1107 1253
## y0[1] 37.23 37.26 4.92 4.60 29.25 45.43 1.00 1907 1810
## y0[2] 19.33 19.40 4.98 4.86 11.36 27.30 1.00 1668 1759
## y0[3] 27.79 27.79 4.83 4.66 20.06 35.60 1.00 1786 1643
## y0[4] 7.35 7.31 5.17 4.85 -0.95 16.11 1.00 784 1357
## y0[5] 25.28 25.30 4.88 4.74 17.55 33.24 1.00 1765 1709
## y0[6] 42.57 42.44 4.98 4.69 34.63 50.81 1.00 1884 1679
## y0[7] 38.74 38.72 5.06 4.78 30.39 46.99 1.00 2097 1950
## y0[8] 22.12 22.08 4.84 4.59 14.11 30.11 1.00 1622 1879
## y0[9] 43.74 43.68 5.19 4.85 35.03 52.44 1.00 2040 2010
## y0[10] 23.51 23.48 4.95 4.51 15.61 31.52 1.00 1702 1872
## y0[11] 4.71 4.56 5.41 5.26 -4.06 13.67 1.00 975 1533
## y0[12] 38.24 38.24 4.86 4.64 30.33 46.45 1.00 1881 1917
## y0[13] 46.05 46.18 5.20 5.03 37.55 54.25 1.00 1897 1798
## y0[14] 32.21 32.19 4.82 4.48 24.26 40.05 1.00 1936 1929
## y0[15] 36.23 36.14 4.94 4.87 28.37 44.33 1.00 1789 1870
## y0[16] 20.39 20.44 4.87 4.77 12.44 28.44 1.00 1782 1929
## y0[17] 24.98 24.96 5.00 4.82 16.95 33.12 1.00 2138 2035
## y0[18] 40.88 40.84 4.99 4.82 32.67 48.97 1.00 2110 2010
## y0[19] 14.22 14.33 4.96 4.68 5.95 22.33 1.00 1364 1541
## y0[20] 31.03 31.12 4.91 4.96 23.09 38.85 1.00 2118 1994
## m1[1] 4.73 4.77 2.61 2.33 0.40 8.92 1.02 273 268
## m1[2] 9.33 9.33 2.21 1.97 5.56 12.82 1.02 272 256
## m1[3] 13.93 13.94 1.84 1.65 10.79 16.80 1.02 276 223
## m1[4] 18.53 18.54 1.52 1.41 15.99 20.94 1.01 294 250
## m1[5] 23.13 23.11 1.26 1.19 21.02 25.15 1.01 358 393
## m1[6] 27.73 27.70 1.12 1.09 25.90 29.59 1.01 592 848
## m1[7] 32.33 32.33 1.15 1.15 30.43 34.19 1.00 1444 1303
## m1[8] 36.93 36.94 1.34 1.36 34.69 39.04 1.00 2092 1520
## m1[9] 41.53 41.54 1.63 1.56 38.79 44.08 1.00 1411 1298
## m1[10] 46.13 46.11 1.98 1.86 42.78 49.26 1.00 967 923
## y1[1] 4.70 4.73 5.39 5.20 -4.18 13.71 1.01 784 1249
## y1[2] 9.23 9.32 5.12 4.90 0.97 17.57 1.00 955 1429
## y1[3] 14.03 14.11 5.09 4.80 5.80 22.60 1.00 1164 1455
## y1[4] 18.49 18.36 5.13 4.90 10.31 26.75 1.00 1459 1828
## y1[5] 23.15 23.26 4.80 4.61 14.95 30.67 1.00 1788 1743
## y1[6] 27.84 27.78 4.90 4.64 20.10 36.04 1.00 1454 1597
## y1[7] 32.30 32.37 4.82 4.61 24.24 40.06 1.00 1833 1971
## y1[8] 36.86 36.81 4.88 4.77 28.67 44.79 1.00 2054 1972
## y1[9] 41.45 41.49 5.00 4.91 33.26 49.55 1.00 2108 2012
## y1[10] 46.16 46.14 5.18 5.00 37.85 54.59 1.00 1872 1976
#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)
mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 2 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 2 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 3 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 4 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 1 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 4 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 3 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 1 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 1 finished in 0.5 seconds.
## Chain 4 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 3 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 3 finished in 0.6 seconds.
## Chain 4 finished in 0.5 seconds.
## Chain 2 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 2 finished in 1.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.7 seconds.
## Total execution time: 1.0 seconds.seeMCMC(mcmc,'[m,x]')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -43.66 -47.24 15.06 11.33 -62.24 -10.23 1.17 17 36
## b0 8.53 8.48 2.18 2.11 4.99 12.16 1.01 722 962
## b1 1.70 1.71 0.15 0.15 1.44 1.95 1.01 676 1234
## s 3.11 3.25 1.63 1.75 0.26 5.59 1.08 52 32
## sx 1.85 1.92 0.95 1.03 0.33 3.34 1.02 88 94
## x0[1] 17.02 16.97 1.16 0.94 15.24 18.92 1.01 1702 1370
## x0[2] 4.09 4.15 2.08 2.52 0.64 6.91 1.02 113 929
## x0[3] 11.41 11.41 1.19 0.90 9.38 13.30 1.02 2083 1154
## x0[4] 2.53 2.43 2.50 3.41 -0.89 6.38 1.06 80 369
## x0[5] 8.43 8.53 1.51 1.68 5.97 10.55 1.01 204 1138
## x0[6] 19.73 19.78 1.24 1.06 17.66 21.72 1.02 1076 1364
## x0[7] 18.38 18.25 1.25 1.09 16.49 20.48 1.01 584 1638
## x0[8] 6.72 6.85 1.52 1.59 4.19 8.85 1.01 241 1054
## x0[9] 21.48 21.28 1.41 1.23 19.49 23.86 1.01 362 1340
## x0[10] 7.97 8.09 1.30 1.19 5.76 9.94 1.01 571 982
## x0[11] -1.71 -1.79 1.39 1.15 -4.00 0.63 1.01 781 1377
## x0[12] 18.10 17.94 1.33 1.19 16.15 20.33 1.01 583 1325
## x0[13] 21.94 21.95 1.30 1.05 19.81 24.08 1.02 1441 1256
## x0[14] 15.71 15.61 1.67 2.07 13.36 18.37 1.02 119 1485
## x0[15] 14.80 14.84 1.55 1.79 12.32 17.06 1.01 198 1374
## x0[16] 6.87 6.92 1.22 0.98 4.78 8.72 1.02 1887 1176
## x0[17] 9.93 9.93 1.19 0.96 7.96 11.78 1.02 1800 1811
## x0[18] 20.47 20.31 1.58 1.68 18.28 23.14 1.01 171 1242
## x0[19] 2.48 2.63 1.45 1.29 0.00 4.55 1.01 556 950
## x0[20] 12.70 12.79 1.19 1.03 10.64 14.58 1.01 1061 1400
## m0[1] 37.45 37.55 1.88 1.51 34.26 40.50 1.04 4144 2460
## m0[2] 15.55 15.52 3.41 4.40 10.90 20.91 1.02 103 1627
## m0[3] 27.93 27.95 1.96 1.49 24.62 31.13 1.05 3573 1781
## m0[4] 12.88 13.18 4.37 5.61 5.66 18.63 1.03 87 725
## m0[5] 22.90 22.85 2.49 2.94 19.38 26.92 1.01 184 1950
## m0[6] 42.04 41.87 2.08 1.74 38.73 45.49 1.02 1288 2140
## m0[7] 39.75 39.97 2.06 1.75 36.27 42.98 1.01 1003 2271
## m0[8] 20.01 19.92 2.49 2.79 16.33 24.02 1.01 212 2618
## m0[9] 45.00 45.22 2.29 2.03 41.03 48.48 1.01 437 2514
## m0[10] 22.13 21.98 2.13 2.03 18.81 25.59 1.01 434 2085
## m0[11] 5.71 5.97 2.41 2.06 1.58 9.49 1.01 944 1799
## m0[12] 39.28 39.41 2.15 1.90 35.71 42.67 1.01 635 2519
## m0[13] 45.78 45.68 2.13 1.85 42.20 49.32 1.04 2592 2314
## m0[14] 35.23 35.33 2.77 3.48 30.80 39.17 1.02 127 2494
## m0[15] 33.70 33.56 2.65 3.12 29.96 38.01 1.01 158 1579
## m0[16] 20.25 20.19 1.97 1.63 17.04 23.42 1.05 3162 2536
## m0[17] 25.44 25.56 1.97 1.57 22.19 28.55 1.03 2369 2122
## m0[18] 43.28 43.43 2.58 2.92 38.82 47.02 1.01 174 2189
## m0[19] 12.82 12.62 2.35 2.16 9.23 16.79 1.01 464 1718
## m0[20] 30.12 29.93 2.00 1.79 26.93 33.41 1.01 1054 2096
## y0[1] 37.45 37.59 3.99 3.10 30.83 44.02 1.03 3997 3174
## y0[2] 15.42 14.53 4.82 4.52 9.02 24.25 1.01 211 2523
## y0[3] 27.82 27.90 4.00 3.02 21.13 34.45 1.03 3706 3046
## y0[4] 12.94 13.98 5.52 5.76 2.86 19.82 1.02 131 1738
## y0[5] 22.79 22.05 4.22 3.37 16.89 30.40 1.01 715 2482
## y0[6] 42.09 41.72 4.15 3.06 35.48 49.21 1.03 2630 1902
## y0[7] 39.76 40.27 4.07 3.17 32.57 45.93 1.02 3946 3147
## y0[8] 20.05 19.25 4.32 3.45 13.83 27.89 1.01 843 2497
## y0[9] 44.98 45.53 4.10 3.26 37.64 51.25 1.01 2032 2933
## y0[10] 22.13 21.58 4.09 3.23 16.10 29.52 1.02 2703 2862
## y0[11] 5.75 6.13 4.20 3.23 -1.46 12.58 1.03 2341 1912
## y0[12] 39.24 39.78 4.16 3.33 31.94 45.63 1.02 2884 2850
## y0[13] 45.77 45.63 4.09 3.18 38.96 52.49 1.04 4118 3323
## y0[14] 35.27 36.09 4.52 3.72 27.09 41.53 1.01 416 3069
## y0[15] 33.81 33.09 4.40 3.65 27.60 41.95 1.01 399 1781
## y0[16] 20.22 20.15 4.01 3.10 13.46 27.13 1.03 3854 2612
## y0[17] 25.47 25.74 4.06 3.23 18.56 32.07 1.03 4008 2482
## y0[18] 43.39 44.16 4.37 3.50 35.31 49.69 1.01 766 3071
## y0[19] 12.71 12.15 4.24 3.30 6.16 20.24 1.01 2123 2426
## y0[20] 30.19 29.79 4.02 3.19 23.81 37.15 1.02 3538 3066
## m1[1] 4.92 4.84 2.46 2.40 0.93 9.05 1.01 703 1036
## m1[2] 9.48 9.44 2.11 2.03 6.07 13.01 1.01 730 1011
## m1[3] 14.04 14.00 1.78 1.70 11.20 16.95 1.01 793 1046
## m1[4] 18.61 18.58 1.48 1.41 16.23 21.05 1.01 922 1164
## m1[5] 23.17 23.16 1.26 1.19 21.13 25.25 1.01 1162 1199
## m1[6] 27.74 27.73 1.15 1.10 25.86 29.64 1.00 1441 1335
## m1[7] 32.30 32.28 1.18 1.17 30.46 34.23 1.00 1499 1231
## m1[8] 36.87 36.86 1.35 1.31 34.75 39.08 1.00 1233 1138
## m1[9] 41.43 41.40 1.60 1.58 38.88 44.06 1.00 1023 1183
## m1[10] 45.99 46.01 1.91 1.89 42.87 49.11 1.00 906 1157
## x1[1] -2.14 -2.12 2.07 1.49 -5.65 1.21 1.01 3796 2739
## x1[2] 0.54 0.55 2.10 1.51 -2.96 4.11 1.01 4025 2613
## x1[3] 3.24 3.25 2.04 1.48 -0.04 6.56 1.01 3691 2511
## x1[4] 6.00 5.95 2.10 1.52 2.66 9.49 1.01 3738 2530
## x1[5] 8.54 8.57 2.07 1.51 5.03 11.94 1.01 4007 3023
## x1[6] 11.28 11.30 2.04 1.52 7.86 14.71 1.01 3872 3034
## x1[7] 13.96 13.96 2.08 1.51 10.50 17.44 1.01 4115 2804
## x1[8] 16.63 16.64 2.05 1.44 13.30 20.09 1.01 3647 2816
## x1[9] 19.39 19.35 2.12 1.47 15.89 23.04 1.01 4045 2954
## x1[10] 22.05 22.04 2.09 1.51 18.64 25.52 1.00 3760 2756
## y1[1] 4.95 4.84 4.27 3.67 -2.07 11.96 1.00 1734 2421
## y1[2] 9.54 9.46 3.99 3.35 3.15 16.48 1.01 1941 2115
## y1[3] 13.97 14.00 3.93 3.26 7.47 20.52 1.01 2543 2807
## y1[4] 18.64 18.64 3.85 3.12 12.24 25.09 1.01 2452 2792
## y1[5] 23.23 23.19 3.70 2.94 17.21 29.42 1.01 3187 3161
## y1[6] 27.72 27.76 3.64 2.85 21.70 33.67 1.01 3655 2929
## y1[7] 32.26 32.27 3.69 2.79 25.99 38.29 1.01 3215 2990
## y1[8] 36.90 36.85 3.75 3.09 30.65 43.24 1.01 3071 2653
## y1[9] 41.45 41.40 3.85 3.19 35.25 47.73 1.01 2953 2574
## y1[10] 46.04 46.02 3.99 3.27 39.44 52.61 1.01 2485 2414
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
outlier
ex10.stan
objective variable have outlier, y~cauchy(b0+b1*x,s)
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~cauchy(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=cauchy_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=cauchy_rng(m1[i],s);
}
}n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -33476
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 26 -32.5138 9.53072e-05 0.000215404 0.8888 0.8888 47
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -32.51
## b0 6.61
## b1 1.82
## s 3.08
## m0[1] 12.08
## m0[2] 15.21
## m0[3] 9.34
## m0[4] 10.17
## m0[5] 18.18
## m0[6] 12.41
## m0[7] 13.50
## m0[8] 19.13
## m0[9] 13.09
## m0[10] 17.96
## m0[11] 16.80
## m0[12] 21.80
## m0[13] 11.14
## m0[14] 16.90
## m0[15] 9.82
## m0[16] 9.75
## m0[17] 18.32
## m0[18] 22.90
## m0[19] 8.19
## m0[20] 8.93
## y0[1] 10.65
## y0[2] 11.30
## y0[3] 12.60
## y0[4] 8.52
## y0[5] 21.04
## y0[6] 9.84
## y0[7] 11.83
## y0[8] 17.60
## y0[9] 14.14
## y0[10] 15.32
## y0[11] 14.40
## y0[12] 19.48
## y0[13] 9.38
## y0[14] 11.98
## y0[15] 14.44
## y0[16] 5.96
## y0[17] 18.74
## y0[18] 24.28
## y0[19] 4.13
## y0[20] 7.25
## m1[1] 8.19
## m1[2] 9.82
## m1[3] 11.46
## m1[4] 13.09
## m1[5] 14.73
## m1[6] 16.36
## m1[7] 17.99
## m1[8] 19.63
## m1[9] 21.26
## m1[10] 22.90
## y1[1] 6.03
## y1[2] 9.41
## y1[3] 15.37
## y1[4] 7.93
## y1[5] 13.85
## y1[6] 12.73
## y1[7] 17.58
## y1[8] 17.76
## y1[9] 24.50
## y1[10] 25.85
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -33.09 -32.73 1.40 1.14 -35.79 -31.58 1.01 439 849
## b0 6.59 6.63 1.70 1.60 3.77 9.28 1.02 287 348
## b1 1.82 1.82 0.35 0.34 1.28 2.41 1.02 349 415
## s 3.52 3.41 0.68 0.60 2.64 4.72 1.00 1493 1305
## m0[1] 12.06 12.08 0.91 0.90 10.54 13.51 1.01 387 552
## m0[2] 15.19 15.22 0.81 0.80 13.82 16.51 1.00 1242 1305
## m0[3] 9.31 9.34 1.26 1.20 7.22 11.26 1.01 299 354
## m0[4] 10.15 10.17 1.14 1.10 8.26 11.92 1.01 310 403
## m0[5] 18.16 18.17 1.08 1.04 16.39 19.93 1.00 1259 1095
## m0[6] 12.39 12.41 0.88 0.88 10.91 13.79 1.01 415 622
## m0[7] 13.47 13.48 0.81 0.79 12.11 14.81 1.00 579 898
## m0[8] 19.11 19.11 1.21 1.19 17.11 21.08 1.01 1063 1101
## m0[9] 13.06 13.08 0.83 0.80 11.66 14.43 1.00 501 749
## m0[10] 17.94 17.94 1.05 1.02 16.21 19.65 1.00 1311 1267
## m0[11] 16.78 16.80 0.92 0.90 15.27 18.26 1.00 1547 1321
## m0[12] 21.78 21.77 1.64 1.61 19.16 24.43 1.01 692 858
## m0[13] 11.12 11.14 1.01 0.99 9.43 12.70 1.01 337 444
## m0[14] 16.87 16.89 0.93 0.90 15.34 18.37 1.00 1536 1321
## m0[15] 9.80 9.82 1.19 1.13 7.80 11.64 1.01 304 382
## m0[16] 9.73 9.75 1.20 1.15 7.72 11.58 1.01 303 382
## m0[17] 18.30 18.32 1.10 1.06 16.50 20.10 1.00 1228 1105
## m0[18] 22.87 22.86 1.83 1.79 19.99 25.83 1.01 624 868
## m0[19] 8.17 8.21 1.44 1.37 5.80 10.42 1.02 289 318
## m0[20] 8.91 8.94 1.32 1.28 6.71 10.98 1.01 295 339
## y0[1] 11.94 11.92 3.69 3.49 5.87 18.10 1.00 1801 1883
## y0[2] 15.14 15.15 3.69 3.50 9.06 20.98 1.00 2109 1819
## y0[3] 9.30 9.48 3.77 3.67 3.24 15.37 1.00 1700 1930
## y0[4] 10.15 10.02 3.75 3.56 4.29 16.53 1.00 1741 2060
## y0[5] 18.13 18.18 3.76 3.55 11.85 24.28 1.00 1907 1750
## y0[6] 12.40 12.43 3.64 3.38 6.35 18.22 1.00 1672 1734
## y0[7] 13.37 13.36 3.74 3.47 7.20 19.72 1.00 1875 1869
## y0[8] 19.14 19.13 3.86 3.79 12.82 25.23 1.00 1816 1907
## y0[9] 13.09 13.11 3.75 3.66 6.93 19.23 1.00 2012 1839
## y0[10] 18.07 18.16 3.72 3.51 11.94 23.95 1.00 1872 1767
## y0[11] 16.82 16.79 3.67 3.48 10.86 22.77 1.00 1822 2023
## y0[12] 21.82 21.69 4.03 3.75 15.40 28.51 1.00 1650 1545
## y0[13] 10.96 10.93 3.73 3.66 4.92 16.95 1.00 1729 1800
## y0[14] 16.88 16.78 3.73 3.58 10.68 23.09 1.00 2014 1683
## y0[15] 9.62 9.69 3.87 3.53 3.49 16.02 1.00 1255 1667
## y0[16] 9.60 9.65 3.85 3.76 3.32 15.80 1.00 1349 1515
## y0[17] 18.30 18.28 3.80 3.68 12.19 24.34 1.00 1808 1931
## y0[18] 23.08 23.15 4.06 3.81 16.39 29.64 1.00 1381 1884
## y0[19] 8.13 8.19 3.88 3.60 1.60 14.59 1.00 1282 1608
## y0[20] 8.99 8.92 3.88 3.77 2.77 15.33 1.00 1264 1872
## m1[1] 8.17 8.21 1.44 1.37 5.80 10.42 1.02 289 318
## m1[2] 9.80 9.83 1.19 1.13 7.80 11.64 1.01 304 382
## m1[3] 11.44 11.46 0.98 0.96 9.80 12.97 1.01 350 469
## m1[4] 13.07 13.08 0.83 0.80 11.67 14.44 1.00 502 749
## m1[5] 14.71 14.73 0.80 0.77 13.38 16.00 1.00 1038 1114
## m1[6] 16.34 16.36 0.88 0.86 14.84 17.75 1.00 1553 1343
## m1[7] 17.97 17.98 1.06 1.03 16.24 19.69 1.00 1302 1220
## m1[8] 19.61 19.59 1.29 1.24 17.51 21.70 1.01 969 1011
## m1[9] 21.24 21.23 1.55 1.53 18.75 23.76 1.01 739 898
## m1[10] 22.87 22.86 1.83 1.79 19.99 25.83 1.01 624 868
## y1[1] 8.13 8.18 3.77 3.52 1.98 14.48 1.00 1209 1768
## y1[2] 9.80 9.82 3.81 3.57 3.35 15.99 1.00 1404 1615
## y1[3] 11.48 11.45 3.68 3.43 5.50 17.72 1.00 1631 1833
## y1[4] 13.04 13.00 3.68 3.54 7.03 19.04 1.00 1570 1796
## y1[5] 14.72 14.65 3.71 3.52 8.60 20.87 1.00 1921 1515
## y1[6] 16.34 16.35 3.67 3.39 10.37 22.43 1.00 1959 1892
## y1[7] 17.99 18.12 3.76 3.46 11.77 24.09 1.00 1930 1812
## y1[8] 19.59 19.54 3.82 3.72 13.38 25.77 1.00 2054 2105
## y1[9] 21.43 21.43 3.89 3.68 15.07 27.86 1.00 1658 1780
## y1[10] 22.89 22.95 4.07 3.91 16.11 29.36 1.00 1536 1588
mdl=cmdstan_model('./ex10.stan')
fn(mdl,data)## Initial log joint probability = -82.5218
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 20 -6.62665 0.000293134 0.000533178 1 1 37
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -6.63
## b0 5.52
## b1 1.89
## s 0.37
## m0[1] 11.20
## m0[2] 14.45
## m0[3] 8.35
## m0[4] 9.22
## m0[5] 17.53
## m0[6] 11.53
## m0[7] 12.66
## m0[8] 18.51
## m0[9] 12.24
## m0[10] 17.29
## m0[11] 16.09
## m0[12] 21.28
## m0[13] 10.22
## m0[14] 16.19
## m0[15] 8.85
## m0[16] 8.78
## m0[17] 17.67
## m0[18] 22.42
## m0[19] 7.16
## m0[20] 7.93
## y0[1] 11.16
## y0[2] 14.58
## y0[3] 6.20
## y0[4] 10.80
## y0[5] 18.49
## y0[6] 13.01
## y0[7] 12.37
## y0[8] 17.91
## y0[9] 12.47
## y0[10] 17.90
## y0[11] 16.49
## y0[12] 21.19
## y0[13] 10.17
## y0[14] 16.21
## y0[15] 8.86
## y0[16] 8.47
## y0[17] 17.78
## y0[18] 21.65
## y0[19] 7.03
## y0[20] 23.14
## m1[1] 7.16
## m1[2] 8.86
## m1[3] 10.55
## m1[4] 12.25
## m1[5] 13.94
## m1[6] 15.64
## m1[7] 17.33
## m1[8] 19.03
## m1[9] 20.72
## m1[10] 22.42
## y1[1] 6.38
## y1[2] 15.73
## y1[3] 2.45
## y1[4] 11.74
## y1[5] 15.35
## y1[6] 14.99
## y1[7] 17.36
## y1[8] 19.29
## y1[9] 20.33
## y1[10] 22.61
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -9.08 -8.74 1.27 1.00 -11.52 -7.73 1.01 458 618
## b0 5.52 5.53 0.29 0.28 5.05 5.98 1.00 696 616
## b1 1.90 1.90 0.07 0.06 1.80 2.01 1.00 821 859
## s 0.45 0.43 0.14 0.13 0.26 0.72 1.01 749 603
## m0[1] 11.21 11.21 0.15 0.15 10.98 11.48 1.00 927 1059
## m0[2] 14.47 14.46 0.16 0.15 14.23 14.75 1.00 1928 1421
## m0[3] 8.36 8.36 0.21 0.21 8.03 8.70 1.00 722 745
## m0[4] 9.23 9.22 0.19 0.19 8.93 9.53 1.00 750 799
## m0[5] 17.56 17.54 0.22 0.20 17.23 17.95 1.00 1697 1536
## m0[6] 11.55 11.55 0.15 0.14 11.32 11.81 1.00 978 1050
## m0[7] 12.68 12.67 0.15 0.14 12.46 12.94 1.00 1278 1365
## m0[8] 18.55 18.53 0.25 0.23 18.18 18.99 1.00 1550 1377
## m0[9] 12.26 12.25 0.15 0.14 12.03 12.51 1.00 1135 1255
## m0[10] 17.32 17.31 0.21 0.20 17.01 17.70 1.00 1737 1488
## m0[11] 16.12 16.10 0.19 0.17 15.84 16.45 1.00 1879 1560
## m0[12] 21.32 21.31 0.33 0.31 20.84 21.91 1.00 1282 1252
## m0[13] 10.23 10.22 0.17 0.16 9.96 10.51 1.00 815 861
## m0[14] 16.22 16.20 0.19 0.18 15.94 16.55 1.00 1869 1542
## m0[15] 8.86 8.86 0.20 0.20 8.55 9.18 1.00 736 784
## m0[16] 8.79 8.79 0.20 0.20 8.47 9.11 1.00 734 782
## m0[17] 17.70 17.69 0.22 0.21 17.37 18.10 1.00 1673 1511
## m0[18] 22.46 22.44 0.36 0.35 21.92 23.10 1.00 1218 1175
## m0[19] 7.17 7.17 0.24 0.24 6.78 7.56 1.00 702 671
## m0[20] 7.94 7.94 0.22 0.22 7.59 8.30 1.00 712 713
## y0[1] 11.38 11.19 23.69 0.65 8.16 14.31 1.00 1883 1828
## y0[2] 14.39 14.45 12.30 0.67 11.65 17.67 1.00 2075 1933
## y0[3] 3.49 8.36 232.58 0.69 5.62 11.42 1.00 1700 1685
## y0[4] 9.28 9.22 12.85 0.71 6.00 12.45 1.00 1900 1728
## y0[5] 17.92 17.58 9.71 0.73 14.95 20.25 1.00 1763 1735
## y0[6] 12.26 11.54 45.29 0.66 7.72 15.07 1.00 2050 1849
## y0[7] 15.00 12.67 146.49 0.69 9.46 15.20 1.00 2092 1833
## y0[8] 18.65 18.54 8.23 0.72 15.60 21.82 1.00 1815 1599
## y0[9] 11.56 12.25 49.36 0.69 9.40 15.56 1.00 1859 1972
## y0[10] 16.92 17.32 52.59 0.68 14.19 20.02 1.00 1920 1950
## y0[11] 15.83 16.10 8.68 0.66 13.05 18.79 1.00 1779 2020
## y0[12] 21.82 21.32 8.55 0.79 18.86 24.49 1.00 1916 1881
## y0[13] 12.31 10.19 115.81 0.70 7.24 13.12 1.00 2002 2041
## y0[14] 16.09 16.21 18.35 0.72 13.01 18.89 1.00 2026 1920
## y0[15] 9.37 8.86 11.58 0.70 5.92 11.60 1.00 1695 1860
## y0[16] 6.87 8.78 72.88 0.70 6.13 11.48 1.00 2018 1857
## y0[17] 18.18 17.71 30.30 0.68 14.76 20.70 1.00 1879 1913
## y0[18] 20.34 22.42 102.03 0.78 19.74 25.42 1.00 2145 1955
## y0[19] 7.08 7.20 7.43 0.69 4.48 9.54 1.00 2062 1943
## y0[20] 8.16 7.94 15.75 0.73 4.77 11.07 1.00 1770 1907
## m1[1] 7.17 7.17 0.24 0.24 6.78 7.56 1.00 702 671
## m1[2] 8.87 8.86 0.20 0.20 8.55 9.18 1.00 736 784
## m1[3] 10.56 10.56 0.16 0.16 10.31 10.84 1.00 844 855
## m1[4] 12.26 12.26 0.15 0.14 12.04 12.52 1.00 1135 1255
## m1[5] 13.96 13.95 0.15 0.15 13.73 14.22 1.00 1790 1399
## m1[6] 15.66 15.65 0.18 0.17 15.39 15.98 1.00 1931 1455
## m1[7] 17.36 17.34 0.22 0.20 17.04 17.74 1.00 1730 1488
## m1[8] 19.06 19.04 0.26 0.24 18.67 19.53 1.00 1489 1365
## m1[9] 20.76 20.75 0.31 0.29 20.30 21.32 1.00 1321 1304
## m1[10] 22.46 22.44 0.36 0.35 21.92 23.10 1.00 1218 1175
## y1[1] 7.25 7.16 23.55 0.69 4.13 10.02 1.00 1618 1669
## y1[2] 8.38 8.83 11.89 0.68 5.50 11.46 1.00 1966 1955
## y1[3] 10.46 10.57 9.25 0.68 7.47 13.45 1.00 1992 1981
## y1[4] 11.99 12.25 33.92 0.67 9.65 14.76 1.00 2023 1853
## y1[5] 14.60 13.95 26.01 0.72 11.26 17.28 1.00 2011 1932
## y1[6] 15.01 15.64 19.49 0.74 12.14 18.33 1.00 1690 1891
## y1[7] 17.08 17.35 24.06 0.67 14.43 20.13 1.00 1824 1710
## y1[8] 19.39 19.03 23.02 0.71 16.06 21.74 1.00 1872 1964
## y1[9] 21.07 20.74 16.73 0.76 18.10 24.02 1.00 1901 1991
## y1[10] 22.34 22.43 12.73 0.79 18.90 25.53 1.00 1791 1745
censored data
objective variable has NA
(x,y) i=1-N
(x0,y0) i=1-N0
x1 i=1-N1, y1=NA
(x,y)~N((mx,my),(sx2,sy2,sxy))
(x0,y0)~N((mx,my),(sx2,sy2,sxy))
x1~N(mx,sx2)ex11-0.stan
data{
int N0;
array[N0] vector[2] xy;
int N1;
vector[N1] x1;
}
parameters{
vector[2] m;
cov_matrix[2] s;
}
model{
target+=multi_normal_lpdf(xy | m, s);
x1~normal(m[1],s[1,1]^.5);
}
generated quantities{
vector[2] xy1;
xy1=multi_normal_rng(m,s);
real cr;
cr=s[1,2]/(s[1,1]*s[2,2])^.5;
}n=30
x=runif(n,0,9)
y=rnorm(n,10+3*x,4)
cor(x,y)## [1] 0.769qplot(x,y)
L=4
n0=sum(x>L)
x0=x[x>L]
y0=y[x>L]
x1=x[x<=L]
n1=sum(x<=L)
cor(x0,y0)## [1] 0.597qplot(x0,y0)
mdl=cmdstan_model('./ex11-0.stan')
data=list(N0=n0,xy=matrix(c(x0,y0),ncol=2),N1=n1,x1=x1)
mle=mdl$optimize(data=data)## Initial log joint probability = -154798
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Exception: multi_normal_lpdf: LDLT_Factor of covariance parameter is not positive definite. last conditional variance is 0. (in '/tmp/Rtmp3zdloC/model-bd3037f1494a.stan', line 12, column 2 to column 39)
## Exception: multi_normal_lpdf: LDLT_Factor of covariance parameter is not positive definite. last conditional variance is 0. (in '/tmp/Rtmp3zdloC/model-bd3037f1494a.stan', line 12, column 2 to column 39)
## 55 -135.468 0.000281555 0.000319153 1 1 98
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -135.47
## m[1] 6.06
## m[2] 27.70
## s[1,1] 4.35
## s[2,1] 10.70
## s[1,2] 10.70
## s[2,2] 50.79
## xy1[1] 6.59
## xy1[2] 34.60
## cr 0.72mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.3 seconds.
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 finished in 0.3 seconds.
## Chain 3 finished in 0.3 seconds.
## Chain 4 finished in 0.3 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.3 seconds.
## Total execution time: 0.4 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -131.26 -130.93 1.71 1.58 -134.62 -129.19 1.00 597 934
## m[1] 6.05 6.04 0.43 0.44 5.35 6.76 1.00 886 938
## m[2] 27.71 27.78 1.66 1.64 24.84 30.26 1.02 457 586
## s[1,1] 5.74 5.39 1.82 1.52 3.53 9.11 1.00 1362 1548
## s[2,1] 14.04 12.84 6.54 5.41 6.02 26.41 1.00 464 631
## s[1,2] 14.04 12.84 6.54 5.41 6.02 26.41 1.00 464 631
## s[2,2] 70.26 63.60 29.40 23.50 36.55 127.07 1.00 502 735
## xy1[1] 6.13 6.12 2.48 2.40 2.11 10.16 1.00 1762 1703
## xy1[2] 27.77 27.80 8.84 8.38 12.90 41.95 1.00 1710 1708
## cr 0.69 0.72 0.14 0.13 0.42 0.87 1.01 460 797
xy=mcmc$draws('xy1')
cor(xy[,,1],xy[,,2])## [1] 0.707qplot(xy[,,1],xy[,,2])
objective variable is censored
y i=1-N, y~N(m,s)
actual ya i=1-Na
lower censored yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
upper censored yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)
cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))
P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)ex11-1.stan
data{
int N;
vector[N] x;
vector[N] y;
real L;
int Nl;
vector[Nl] xl;
real U;
int Nu;
vector[Nu] xu;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
for(i in 1:Nl)
target+=normal_lcdf(L | b0+b1*xl[i],s);
for(i in 1:Nu)
target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)
xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcitionmdl=cmdstan_model('./ex8-0.stan')
data=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)
mle=mdl$optimize(data=data)## Initial log joint probability = -3560.18
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 20 -23.3087 0.00437841 0.000468356 1 1 37
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -23.31
## b0 17.68
## b1 1.48
## s 3.64
## m0[1] 24.66
## m0[2] 21.36
## m0[3] 23.03
## m0[4] 26.06
## m0[5] 25.95
## m0[6] 20.02
## m0[7] 21.26
## m0[8] 22.64
## m0[9] 21.12
## m0[10] 25.26
## m0[11] 20.54
## m0[12] 24.27
## m0[13] 28.13
## y0[1] 18.42
## y0[2] 13.38
## y0[3] 20.01
## y0[4] 26.45
## y0[5] 25.72
## y0[6] 19.16
## y0[7] 25.18
## y0[8] 27.03
## y0[9] 25.98
## y0[10] 22.96
## y0[11] 18.84
## y0[12] 24.86
## y0[13] 29.96
## m1[1] 17.86
## m1[2] 19.25
## m1[3] 20.65
## m1[4] 22.04
## m1[5] 23.44
## m1[6] 24.84
## m1[7] 26.23
## m1[8] 27.63
## m1[9] 29.02
## m1[10] 30.42
## y1[1] 13.67
## y1[2] 21.14
## y1[3] 21.89
## y1[4] 28.13
## y1[5] 28.91
## y1[6] 29.89
## y1[7] 25.31
## y1[8] 28.42
## y1[9] 23.86
## y1[10] 30.50mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -23.64 -23.29 1.31 1.07 -26.18 -22.18 1.00 585 809
## b0 17.51 17.35 3.15 3.09 12.76 22.97 1.01 309 285
## b1 1.53 1.56 0.74 0.74 0.25 2.68 1.00 361 494
## s 4.46 4.29 1.10 0.96 3.06 6.38 1.00 1402 1344
## m0[1] 24.72 24.69 1.39 1.32 22.52 27.03 1.00 1797 1326
## m0[2] 21.32 21.24 1.62 1.51 18.79 24.02 1.01 378 516
## m0[3] 23.03 23.02 1.26 1.14 20.97 25.06 1.00 780 929
## m0[4] 26.17 26.19 1.82 1.67 23.27 29.10 1.00 1097 1151
## m0[5] 26.05 26.07 1.78 1.64 23.19 28.92 1.00 1144 1137
## m0[6] 19.93 19.81 2.12 2.00 16.71 23.44 1.01 324 333
## m0[7] 21.21 21.14 1.66 1.54 18.64 23.98 1.01 370 486
## m0[8] 22.63 22.59 1.31 1.18 20.49 24.76 1.01 605 786
## m0[9] 21.06 20.97 1.70 1.59 18.41 23.94 1.01 361 491
## m0[10] 25.34 25.36 1.55 1.44 22.86 27.90 1.00 1540 1199
## m0[11] 20.46 20.34 1.91 1.78 17.52 23.68 1.01 337 364
## m0[12] 24.32 24.32 1.32 1.23 22.20 26.49 1.00 1822 1367
## m0[13] 28.30 28.33 2.67 2.57 23.84 32.56 1.00 672 931
## y0[1] 24.75 24.76 4.82 4.64 16.87 32.84 1.00 2079 1973
## y0[2] 21.16 21.12 4.85 4.52 13.48 28.91 1.00 1539 1578
## y0[3] 22.75 22.81 4.68 4.41 15.17 30.32 1.00 1961 1910
## y0[4] 26.18 26.28 4.99 4.75 17.96 34.30 1.00 1690 1822
## y0[5] 26.10 26.11 4.79 4.48 18.73 33.92 1.00 1966 1822
## y0[6] 20.06 19.93 5.01 4.83 12.13 28.07 1.00 864 1370
## y0[7] 21.10 21.10 4.90 4.57 13.37 28.95 1.00 1248 1766
## y0[8] 22.48 22.37 4.66 4.37 14.91 29.90 1.00 1967 1971
## y0[9] 21.03 21.12 4.89 4.52 13.16 28.96 1.00 1030 1753
## y0[10] 25.33 25.15 4.87 4.47 17.59 33.66 1.00 1979 1846
## y0[11] 20.38 20.21 4.94 4.59 12.54 28.50 1.00 1060 1492
## y0[12] 24.21 24.17 4.99 4.61 16.30 32.43 1.00 1953 1884
## y0[13] 28.28 28.50 5.33 4.91 19.29 36.87 1.00 1343 1538
## m1[1] 17.69 17.55 3.07 3.01 13.07 23.05 1.01 309 278
## m1[2] 19.13 18.99 2.44 2.37 15.41 23.18 1.01 315 290
## m1[3] 20.58 20.46 1.87 1.75 17.68 23.74 1.01 340 375
## m1[4] 22.02 21.96 1.43 1.28 19.73 24.35 1.01 459 681
## m1[5] 23.46 23.45 1.25 1.13 21.40 25.55 1.00 1102 1270
## m1[6] 24.90 24.90 1.43 1.36 22.63 27.28 1.00 1766 1326
## m1[7] 26.34 26.38 1.88 1.74 23.34 29.34 1.00 1033 1165
## m1[8] 27.78 27.81 2.45 2.34 23.73 31.71 1.00 727 958
## m1[9] 29.23 29.27 3.08 2.91 24.06 34.00 1.00 605 868
## m1[10] 30.67 30.69 3.73 3.56 24.37 36.53 1.00 531 740
## y1[1] 17.66 17.56 5.53 5.12 8.79 26.98 1.01 826 1236
## y1[2] 19.18 19.11 5.32 5.02 10.69 27.95 1.00 960 1510
## y1[3] 20.73 20.72 5.03 4.51 12.68 29.07 1.00 999 1256
## y1[4] 22.06 22.04 4.99 4.49 14.23 30.00 1.00 1480 1974
## y1[5] 23.46 23.47 4.61 4.39 16.30 30.75 1.00 2163 2037
## y1[6] 25.01 24.92 4.78 4.37 17.23 32.81 1.00 2107 1934
## y1[7] 26.27 26.42 5.07 4.79 18.12 34.59 1.00 1601 1777
## y1[8] 27.68 27.61 5.32 5.05 19.18 36.34 1.00 1430 1856
## y1[9] 29.21 29.28 5.50 5.13 20.39 37.90 1.00 1213 1407
## y1[10] 30.70 30.77 5.95 5.59 20.97 40.36 1.00 879 1461
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
data=list(N=n,x=xya$x,y=xya$y,
L=L,Nl=nl,xl=xyl$x,
U=U,Nu=nu,xu=xyu$x,
N1=n1,x1=x1)
mdl=cmdstan_model('./ex11-1.stan')
mle=mdl$optimize(data=data)## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Initial log joint probability = -117.14
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 31 -29.0816 0.000876035 4.91276e-07 1 1 41
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -29.08
## b0 12.65
## b1 2.78
## s 4.30
## m0[1] 25.78
## m0[2] 19.58
## m0[3] 22.71
## m0[4] 28.41
## m0[5] 28.20
## m0[6] 17.06
## m0[7] 19.39
## m0[8] 21.98
## m0[9] 19.12
## m0[10] 26.91
## m0[11] 18.03
## m0[12] 25.05
## m0[13] 32.30
## y0[1] 25.50
## y0[2] 26.64
## y0[3] 25.35
## y0[4] 31.36
## y0[5] 24.37
## y0[6] 20.58
## y0[7] 15.76
## y0[8] 21.35
## y0[9] 16.90
## y0[10] 29.00
## y0[11] 17.16
## y0[12] 30.25
## y0[13] 34.52
## m1[1] 12.99
## m1[2] 15.61
## m1[3] 18.24
## m1[4] 20.86
## m1[5] 23.48
## m1[6] 26.11
## m1[7] 28.73
## m1[8] 31.35
## m1[9] 33.98
## m1[10] 36.60
## y1[1] 12.40
## y1[2] 9.25
## y1[3] 23.55
## y1[4] 22.33
## y1[5] 30.85
## y1[6] 33.42
## y1[7] 19.85
## y1[8] 37.26
## y1[9] 25.63
## y1[10] 39.91mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc,'m',ch=T)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -29.25 -28.90 1.36 1.08 -31.80 -27.79 1.01 564 616
## b0 11.71 11.91 3.13 3.01 6.40 16.34 1.02 241 322
## b1 3.04 2.98 0.71 0.68 2.02 4.25 1.01 277 437
## s 5.33 5.04 1.40 1.15 3.63 8.03 1.00 1099 1073
## m0[1] 26.07 25.98 1.41 1.32 23.91 28.51 1.00 1554 1214
## m0[2] 19.29 19.38 1.71 1.56 16.51 21.95 1.01 319 402
## m0[3] 22.71 22.69 1.35 1.25 20.43 24.90 1.01 586 803
## m0[4] 28.94 28.82 1.76 1.63 26.28 32.06 1.00 1350 1029
## m0[5] 28.72 28.60 1.72 1.59 26.11 31.76 1.00 1394 1120
## m0[6] 16.53 16.63 2.17 2.04 12.93 19.81 1.02 262 351
## m0[7] 19.08 19.16 1.74 1.60 16.22 21.78 1.01 312 400
## m0[8] 21.91 21.92 1.40 1.27 19.58 24.17 1.01 477 717
## m0[9] 18.78 18.87 1.78 1.62 15.86 21.51 1.01 303 381
## m0[10] 27.30 27.20 1.53 1.42 24.99 30.02 1.00 1624 1233
## m0[11] 17.59 17.70 1.98 1.87 14.31 20.59 1.01 276 364
## m0[12] 25.27 25.20 1.36 1.26 23.18 27.47 1.00 1339 1262
## m0[13] 33.20 33.02 2.52 2.31 29.52 37.93 1.00 624 970
## y0[1] 25.86 25.78 5.84 5.42 16.32 35.41 1.00 2094 1931
## y0[2] 19.17 19.26 5.73 5.32 9.81 28.21 1.00 1464 1673
## y0[3] 22.59 22.61 5.68 5.42 13.44 31.64 1.00 1894 1856
## y0[4] 29.18 28.96 5.96 5.51 19.57 38.99 1.00 1753 1962
## y0[5] 28.78 28.86 5.85 5.57 19.90 38.24 1.00 1818 1921
## y0[6] 16.63 16.73 6.12 5.43 6.71 26.64 1.00 1357 1450
## y0[7] 19.33 19.24 5.87 5.22 10.02 28.80 1.01 1261 1528
## y0[8] 22.02 21.96 5.70 5.24 12.90 31.12 1.00 1746 1664
## y0[9] 18.67 18.59 5.92 5.42 9.67 28.19 1.00 1559 1601
## y0[10] 27.37 27.24 5.56 5.19 18.51 36.55 1.00 1758 1860
## y0[11] 17.58 17.75 5.95 5.26 7.49 26.91 1.00 1459 1560
## y0[12] 25.17 25.29 5.55 5.10 16.27 34.38 1.00 2022 2014
## y0[13] 33.40 33.36 6.08 5.90 23.70 43.19 1.00 1287 1678
## m1[1] 12.08 12.28 3.05 2.93 6.89 16.59 1.02 241 321
## m1[2] 14.95 15.10 2.47 2.36 10.84 18.61 1.02 250 319
## m1[3] 17.82 17.91 1.94 1.83 14.59 20.77 1.01 281 364
## m1[4] 20.69 20.74 1.52 1.37 18.23 23.12 1.01 379 518
## m1[5] 23.56 23.51 1.32 1.22 21.41 25.69 1.00 760 1087
## m1[6] 26.43 26.34 1.44 1.35 24.23 28.95 1.00 1617 1235
## m1[7] 29.30 29.20 1.81 1.68 26.56 32.50 1.00 1288 1039
## m1[8] 32.17 32.02 2.32 2.12 28.78 36.48 1.00 726 1054
## m1[9] 35.04 34.82 2.89 2.69 30.83 40.38 1.00 521 940
## m1[10] 37.90 37.62 3.50 3.29 32.78 44.27 1.00 440 791
## y1[1] 12.01 12.25 6.23 5.83 1.40 21.80 1.00 707 1196
## y1[2] 14.91 15.39 6.06 5.60 4.34 23.96 1.00 989 1122
## y1[3] 17.84 17.85 5.69 5.27 8.50 27.11 1.00 1458 1842
## y1[4] 20.69 20.70 5.56 5.09 11.59 29.64 1.00 1884 1878
## y1[5] 23.53 23.59 5.66 5.14 14.23 32.81 1.00 1870 1958
## y1[6] 26.42 26.40 5.71 5.20 17.35 35.91 1.00 1954 1926
## y1[7] 29.15 29.15 5.93 5.39 19.40 39.06 1.00 1952 1957
## y1[8] 32.13 31.89 5.89 5.40 22.96 41.91 1.00 1656 1848
## y1[9] 34.95 34.72 6.25 5.79 25.34 45.69 1.00 1652 1561
## y1[10] 37.84 37.51 6.70 6.29 27.67 48.79 1.00 1169 1487
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
sensitivity/specificity analysis
ex14.stan
estimating sens and spec
data {
int N;
array[N] int x;
array[N] int y;
}
parameters {
real<lower=0,upper=1> p;
real<lower=0,upper=1> se;
real<lower=0,upper=1> sp;
}
model {
p~uniform(0,1);
se~uniform(0,1);
sp~uniform(0,1);
for (i in 1:N) {
y[i]~bernoulli(x[i]*se+(1-x[i])*(1-sp));
}
}
generated quantities {
real ppv;
real npv;
ppv=se*p/((se*p)+((1-p)*(1-sp)));
npv=(1-p)*sp/(((1-p)*sp)+(p*(1-se)));
}n=20
data=list(N=n,
x=sample(0:1,n,replace=T),
y=sample(0:1,n,replace=T))
mdl=cmdstan_model('./ex14.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -13.0211
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 5 -12.0405 0.000152385 1.84742e-05 0.9347 0.9347 8
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -12.04
## p 0.70
## se 0.75
## sp 0.37
## ppv 0.73
## npv 0.40mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -18.13 -17.81 1.32 1.12 -20.67 -16.70 1.00 708 896
## p 0.49 0.48 0.29 0.37 0.04 0.95 1.00 673 642
## se 0.71 0.73 0.12 0.12 0.50 0.89 1.00 2566 1340
## sp 0.41 0.40 0.14 0.15 0.18 0.65 1.00 2449 1411
## ppv 0.52 0.53 0.29 0.37 0.05 0.95 1.01 685 667
## npv 0.57 0.62 0.29 0.36 0.07 0.97 1.00 734 579
2x2 cross table
Effect occur y=1 with probabilty p01, p11 from each cause x{0,1}
event frequncy nxy of effect y{0,1} by cause x{0,1}
n01~B(n0.,p0)
n11~B(n1.,p1)
n01=n0p0, n00=n0(1-p0)
n11=n1p1, n10=n1(1-p1)
p00=n00/n=n0(1-p0)/n, p01=n01/n=n0p0/n
p10=n10/n=n1(1-p1)/n, p11=n11/n=n1p1/n
Cramer'V (chi2/n/(min(row,column)-1))^.5
in 2x2
crv =(n11*n00-n10*n01)/(n0.*n1.*n.0*n.1)^.5
=(n0(1-p0)n1p1-n0p0n1(1-p1))/(n0n1(n0(1-p0)+n1(1-p1))(n0p0n1p1))^.5
kappa coefficient k=(po-pe)/(1-pe)
po: Observed agreement (proportion of times both raters agreed)
pe: Expected agreement under independence
po=p00+p11
=(n0(1-p0)+n1p1)/n
pe=(p00+p01)(p00+p10)(p10+p11)(p01+p11)
=n0/n*(n0(1-p0)+n1(1-p1))/n*(n0p0+n1p1)/n*n1/nex16-1.stan
data {
int n;
int n0;
int n01;
int n1;
int n11;
}
parameters {
real<lower=0,upper=1> p0;
real<lower=0,upper=1> p1;
}
model {
n01~binomial(n0,p0);
n11~binomial(n1,p1);
}
generated quantities {
real RR;
RR=p1/p0;
real OR;
OR=(p1/(1-p1))/(p0/(1-p0));
}ex16-2.stan
data {
int n;
int n0;
int n01;
int n1;
int n11;
}
parameters {
real<lower=0,upper=1> p0;
real<lower=0,upper=1> p1;
}
model {
n01~binomial(n0,p0);
n11~binomial(n1,p1);
}
generated quantities {
real RR;
RR=p1/p0;
real OR;
OR=(p1/(1-p1))/(p0/(1-p0));
real crv;
crv=(n0*(1-p0)*n1*p1-n0*p0*n1*(1-p1))/(n0*n1*(n0*(1-p0)+n1*(1-p1))*(n0*p0+n1*p1))^.5;
real k;
real po;
po=(n0*(1-p0)+n1*p1)/n;
real pe;
pe=n0/n*(n0*(1-p0)+n1*(1-p1))/n*(n0*p0+n1*p1)/n*n1/n;
k=(po-pe)/(1-pe);
}n0=30
n01=rbinom(1,n0,0.3)
n1=30
n11=rbinom(1,n1,0.6)
data=list(n=n0+n1,n0=n0,n01=n01,n1=n1,n11=n11)
mdl=cmdstan_model('./ex16-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -71.8484
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 5 -39.2858 0.00286964 0.000356069 0.9698 0.9698 8
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -39.29
## p0 0.33
## p1 0.60
## RR 1.80
## OR 3.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -43.24 -42.93 1.04 0.72 -45.32 -42.27 1.00 942 1027
## p0 0.34 0.34 0.08 0.09 0.21 0.48 1.00 1966 1351
## p1 0.59 0.59 0.08 0.09 0.46 0.73 1.00 2150 1342
## RR 1.85 1.75 0.60 0.52 1.11 2.91 1.00 1986 1453
## OR 3.34 2.88 1.98 1.55 1.23 6.95 1.00 1976 1473
data=list(n=n0+n1,n0=n0,n01=n01,n1=n1,n11=n11)
mdl=cmdstan_model('./ex16-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -47.5796
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 4 -39.2858 0.000782113 1.21406e-05 1 1 7
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -39.29
## p0 0.33
## p1 0.60
## RR 1.80
## OR 3.00
## crv 0.27
## k 0.61
## po 0.63
## pe 0.06mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -43.24 -42.93 1.04 0.72 -45.32 -42.27 1.00 942 1027
## p0 0.34 0.34 0.08 0.09 0.21 0.48 1.00 1966 1351
## p1 0.59 0.59 0.08 0.09 0.46 0.73 1.00 2150 1342
## RR 1.85 1.75 0.60 0.52 1.11 2.91 1.00 1986 1453
## OR 3.34 2.88 1.98 1.55 1.23 6.95 1.00 1976 1473
## crv 0.25 0.26 0.12 0.13 0.05 0.44 1.00 1976 1486
## k 0.60 0.60 0.06 0.07 0.49 0.70 1.00 1975 1473
## po 0.62 0.63 0.06 0.06 0.53 0.72 1.00 1976 1450
## pe 0.06 0.06 0.00 0.00 0.06 0.06 1.00 1519 1277
bimodal distribution, mixed normal distribution
n=100
y0=rnorm(n,0,1)
y1=rnorm(n,-5,1)
y2=rnorm(n,5,1)
y=sample(c(y0,y1,y2),n)
density(y) |> plot()
EM algorithm
library(mclust)
rst=Mclust(y)
summary(rst)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust E (univariate, equal variance) model with 4 components:
##
## log-likelihood n df BIC ICL
## -230 100 8 -497 -500
##
## Clustering table:
## 1 2 3 4
## 36 24 7 33rst$parameters## $pro
## [1] 0.3600 0.2415 0.0664 0.3321
##
## $mean
## 1 2 3 4
## -5.32 -0.14 3.09 5.44
##
## $variance
## $variance$modelName
## [1] "E"
##
## $variance$d
## [1] 1
##
## $variance$G
## [1] 4
##
## $variance$sigmasq
## [1] 0.51plot(rst)
ex17-1.stan
data {
int N;
vector[N] y;;
}
parameters {
simplex[3] h; //ratio of mix
real m1;
real m2;
real m3;
real<lower=0> s1;
real<lower=0> s2;
real<lower=0> s3;
}
model {
s1~cauchy(0,5);
s2~cauchy(0,5);
s3~cauchy(0,5);
for (i in 1:N) {
vector[3] p;
p[1]=log(h[1]) + normal_lpdf(y[i] | m1, s1);
p[2]=log(h[2]) + normal_lpdf(y[i] | m2, s2);
p[3]=log(h[3]) + normal_lpdf(y[i] | m3, s3);
target+=log_sum_exp(p);
}
}mdl=cmdstan_model('./ex17-1.stan')
data=list(N=n,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -474.65
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 93 -289.352 0.000914839 0.00664515 0.8877 0.8877 137
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -289.35
## h[1] 0.92
## h[2] 0.05
## h[3] 0.03
## m1 0.08
## m2 -0.38
## m3 0.31
## s1 4.77
## s2 0.08
## s3 0.03mcmc=goMCMC(mdl,data,smp=2000)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 2100 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 2100 [ 4%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 2100 [ 0%] (Warmup)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 2100 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 2100 [ 4%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 2100 [ 0%] (Warmup)
## Chain 1 Iteration: 1100 / 2100 [ 52%] (Sampling)
## Chain 3 Iteration: 1100 / 2100 [ 52%] (Sampling)
## Chain 2 Iteration: 101 / 2100 [ 4%] (Sampling)
## Chain 1 Iteration: 2100 / 2100 [100%] (Sampling)
## Chain 4 Iteration: 101 / 2100 [ 4%] (Sampling)
## Chain 1 finished in 1.3 seconds.
## Chain 3 Iteration: 2100 / 2100 [100%] (Sampling)
## Chain 3 finished in 1.4 seconds.
## Chain 4 Iteration: 1100 / 2100 [ 52%] (Sampling)
## Chain 2 Iteration: 1100 / 2100 [ 52%] (Sampling)
## Chain 4 Iteration: 2100 / 2100 [100%] (Sampling)
## Chain 4 finished in 28.3 seconds.
## Chain 2 Iteration: 2100 / 2100 [100%] (Sampling)
## Chain 2 finished in 31.6 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 15.7 seconds.
## Total execution time: 31.8 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -257.32 -260.91 15.38 23.12 -276.31 -239.76 2.02 5 35
## h[1] 0.13 0.10 0.13 0.14 0.00 0.32 1.73 6 128
## h[2] 0.35 0.35 0.06 0.06 0.26 0.44 1.14 18 56
## h[3] 0.51 0.52 0.15 0.21 0.31 0.73 1.81 5 43
## m1 614.14 1.77 863.80 4.51 -0.30 2522.67 2.17 5 13
## m2 -0.04 -0.27 5.29 7.76 -5.50 5.50 1.90 5 36
## m3 0.10 0.50 4.17 6.10 -5.43 5.30 2.84 4 26
## s1 16.09 1.10 388.25 0.71 0.58 30.95 1.41 10 159
## s2 0.84 0.81 0.21 0.17 0.51 1.22 1.43 8 28
## s3 2.06 1.93 1.18 1.59 0.70 3.83 2.16 5 30
